Preamble

# BASH Interface packages
import os
import glob
import subprocess

# Dataframe pacakegs
import pandas as pd
import numpy as np

# Jupyter plugins
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import warnings
import math

# Plotting packages
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
warnings.filterwarnings("ignore")

# Theme stuff
# from jupyterthemes import jtplot
# jtplot.style(theme='grade3')
# jtplot.style(context='talk', fscale=1.4, spines=False, gridlines='--')
# jtplot.style(ticks=True, grid=False, figsize=(6, 4.5))

Processing

Def functions and variables

Base

# Defines a few functions that:

# Filters through the pd dataset
def filter_class(basename, ionname, energy):
    return(list(data[data["base_name"] == basename][data["ion_class"] == ionname][energy]))

def filter_class2(basename, ionname, energy):
    return(list(homodata[homodata["base_name"] == basename][homodata["ion_class"] == ionname][energy]))

# Runs bash commands 
def runbash(command):
    proc = subprocess.run(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, check=False, text=True, executable='/bin/bash')
    output = proc.stdout.strip('\n')
    if output != "":
        return(output)
    
def Eh2eV(df, column):
    for icolumn in column:
        for i in range(len(df[icolumn])):
            df[icolumn][i] = df[icolumn][i] * 27.2116

Pairs

def filter_class_pairs(basename, ionname, energy):
    return(list(datapairs[datapairs["base_name"] == basename][datapairs["ion_name"] == ionname][energy]))

def filter_class2_pairs(basename, ionname, energy):
    return(list(homodatapairs[homodatapairs["base_name"] == basename][homodatapairs["ion_name"] == ionname][energy]))

def filter_ions_pairs(ionname, energy):
    return(list(datapairs[datapairs["ion_name"] == ionname][energy]))

Globals

# The lists to cycle through
basenames = ["ethane", "ethene", "ethyne", "benzene", "furan", "pyrrole"]
ionnames = ["BF2(CF3)2", "BF4", "C(CN)3", "N(SOOCF3)2", "Cl" ,"N(NC)2", "C1MPyr", "ImidMe2", "NMe4", "PMe4", "Na", "Li", "Pd", "Ni", "tempo"]
classnames = ["cation", "anion", "metal"]
classnamespairs = ["high", "low"]
energynames = ["electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "total", 'dhomo', 'dlumo', 'homo', 'lumo', 'gap', 'nbo_homo', 'nbo_lumo', 'nbo_gap', 'dplusel']

(legacy) Import and process data

Since removing from the server it was running on with complptete access to live calcualtions as they're coming in, the data has been exported to csv files to be imported cleanly.

I could have converted these blocks to raw, but then my code hiding extension wouldn't hide them.

Base

# homodata = pd.read_csv("./IonsData/homo-lumo-base.csv")


# # Sets the dataframe column type to float
# homodata.homo = homodata.homo.astype(float)
# homodata.lumo = homodata.lumo.astype(float)

# lone = pd.DataFrame(columns=['base_name', 'ion_class', 'ion_name', 'homo', 'lumo'])
# lone = lone.append(homodata[homodata["base_name"] == "alkyl"], ignore_index = True)
# lone = lone.append(homodata[homodata["base_name"] == "anions"], ignore_index = True)
# lone = lone.append(homodata[homodata["base_name"] == "cations"], ignore_index = True)
# lone = lone.append(homodata[homodata["base_name"] == "aromatic"], ignore_index = True)

# droplist = []
# for i in reversed(range(len(list(homodata.base_name)))):
#     if list(homodata.base_name)[i] == "alkyl" or list(homodata.base_name)[i] == "aromatic" or list(homodata.base_name)[i] == "cations" or list(homodata.base_name)[i] == "anions":
#         droplist = droplist + [i]
        
# homodata = homodata.drop(droplist)

# # Initalise the dataframe
# data = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'electrostatics', 'exchange', 'induction', 'dispersion', 'total', 'charge_transfer', 'dhomo', 'dlumo', 'bhomo', 'blumo', 'plotsize'])

# # This runs through the files and pulls out important data
# infiles = glob.glob("../4-SAPT/*/*.out")
# for i in infiles:
#     if "slurm" not in i:
        
#         fullfilepath = runbash("readlink -f '" + i + "'")
#         filename = runbash("/bin/basename -s .out '" + fullfilepath + "'")
#         ion_name = runbash("echo '" + filename + "' | cut -d'-' -f2")
#         base_name = runbash("echo '" + filename + "' | cut -d'-' -f1")
        
#         sort = {"ethane"  :1, "ethene" : 2, "ethyne" : 3, "benzene" : 4, "furan" : 5, "pyrrole" : 6}
        
#         classes = [["BF2(CF3)2", "BF4", "C(CN)3", "N(SOOCF3)2", "Cl" ,"N(NC)2", "tempo"],["C1MPyr", "ImidMe2", "NMe4", "PMe4", "Na", "Li"],["Pd", "Ni"]]
#         if ion_name in classes[0]:
#             ion_class = "anion"
#         elif ion_name in classes[1]:
#             ion_class = "cation"
#         elif ion_name in classes[2]:
#             ion_class = "metal"
    
#         electrostatics = runbash("cat '" + i + "' | grep 'Electrostatics  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         exchange  = runbash("cat '" + i + "' | grep 'Exchange  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         induction  = runbash("cat '" + i + "' | grep 'Induction  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         dispersion  = runbash("cat '" + i + "' | grep 'Dispersion  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         total  = runbash("cat '" + i + "' | grep ' Total SAPT2+  ' -m1 | tr -s ' ' | cut -d' ' -f8")
#         charge_transfer  = runbash("cat '" + i + "' | grep 'SAPT Charge Transfer' | tail -n 1 | tr -s ' ' | cut -d' ' -f9")
        
#         newdata = {"sort" : sort[base_name], "base_name" : base_name, "ion_class" : ion_class, "ion_name" : ion_name, "electrostatics" : electrostatics, "exchange" : exchange, "induction" : induction, "dispersion" : dispersion, "total" : total, "charge_transfer" : charge_transfer, "dhomo" : 0, "dlumo" : 0, "plotsize" : 0}
#         data = data.append(newdata, ignore_index=True)

# # Drops any frames that have a "None" value in charge_trasfer (the last SAPT calculation)

# data = data[data['charge_transfer'].isnull() == False]

# # Sets the dataframe column type to float
# data.electrostatics = data.electrostatics.astype(float)
# data.exchange = data.exchange.astype(float)
# data.induction = data.induction.astype(float)
# data.dispersion = data.dispersion.astype(float)
# data.total = data.total.astype(float)
# data.charge_transfer = data.charge_transfer.astype(float)
# data.dhomo = data.dhomo.astype(float)
# data.dlumo = data.dlumo.astype(float)
# data.bhomo = data.bhomo.astype(float)
# data.blumo = data.blumo.astype(float)
# data.plotsize = data.plotsize.astype(int)

# data = pd.merge(data, homodata, on=['base_name', 'ion_name', 'ion_class'])

# warnings.filterwarnings("ignore")
# for i in range(len(data.base_name)):
#     for j in range(len(lone.homo)):

#         if data.base_name[i] == lone.ion_name[j]:
#             data.dhomo[i] = data.homo[i] - lone.homo[j]
#             data.dlumo[i] = data.lumo[i] - lone.lumo[j]
#             data.bhomo[i] = lone.homo[j]
#             data.blumo[i] = lone.lumo[j]

# for i in range(len(data.base_name)):
#     if math.isnan(data.total[i]) == False:
#         data.plotsize[i] = abs(int(data.total[i]))
    
# Eh2eV(data, ["lumo", "homo", "dhomo", "dlumo", "bhomo", "blumo"])
# Eh2eV(lone, ["lumo", "homo"])

Pairs

# lonedata = pd.read_csv("./IonsData/homo-lumo-base.csv")

# lonepairs = pd.DataFrame(columns=['base_name', 'ion_class', 'ion_name', 'homo', 'lumo', 'gap', 'nbo_homo', 'nbo_lumo', 'nbo_gap'])
# lonepairs = lonepairs.append(lonedata[lonedata["base_name"] == "alkyl"], ignore_index = True)
# lonepairs = lonepairs.append(lonedata[lonedata["base_name"] == "aromatic"], ignore_index = True)


# # Sets the dataframe column type to float
# lonepairs.homo = lonepairs.homo.astype(float)
# lonepairs.lumo = lonepairs.lumo.astype(float)

# lonepairs.drop('ion_class', axis=1, inplace=True)

# for i in range(len(lonepairs)):
#     lonepairs.base_name[i] = lonepairs.ion_name[i]
#     lonepairs.ion_name[i] = "lone"


# homodatapairs = pd.read_csv("./IonsData/homo-lumo-pairs.csv")

# # Sets the dataframe column type to float
# homodatapairs.homo = homodatapairs.homo.astype(float)
# homodatapairs.lumo = homodatapairs.lumo.astype(float)

# homodatapairs.drop('ion_class', axis=1, inplace=True)

# # Initalise the dataframe
# datapairs = pd.DataFrame(columns=['sort', 'base_name', 'ion_name', 'electrostatics', 'exchange', 'induction', 'dispersion', 'total', 'charge_transfer', 'dhomo', 'dlumo', 'bhomo', 'blumo', 'plotsize'])


# # This runs through the files and pulls out important data
# infiles = glob.glob("../7-IonPairsSAPT/*.out")
# for i in infiles:
#     if "slurm" not in i:
        
#         fullfilepath = runbash("readlink -f '" + i + "'")
#         filename = runbash("/bin/basename -s .out '" + fullfilepath + "'")
#         ion_name = runbash("echo '" + filename + "' | cut -d'-' -f2")
#         base_name = runbash("echo '" + filename + "' | cut -d'-' -f1")

#         electrostatics = runbash("cat '" + i + "' | grep 'Electrostatics  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         exchange  = runbash("cat '" + i + "' | grep 'Exchange  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         induction  = runbash("cat '" + i + "' | grep 'Induction  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         dispersion  = runbash("cat '" + i + "' | grep 'Dispersion  ' -m1 | tr -s ' ' | cut -d' ' -f7")
#         total  = runbash("cat '" + i + "' | grep ' Total SAPT2+  ' -m1 | tr -s ' ' | cut -d' ' -f8")
#         charge_transfer  = runbash("cat '" + i + "' | grep 'SAPT Charge Transfer' | tail -n 1 | tr -s ' ' | cut -d' ' -f9")
        
#         newdata = {"sort" : sort[base_name], "base_name" : base_name, "ion_name" : ion_name, "electrostatics" : electrostatics, "exchange" : exchange, "induction" : induction, "dispersion" : dispersion, "total" : total, "charge_transfer" : charge_transfer, "dhomo" : 0, "dlumo" : 0, "plotsize" : 0}
#         datapairs = datapairs.append(newdata, ignore_index=True)

# # Drops any frames that have a "None" value in charge_trasfer (the last SAPT calculation)

# datapairs = datapairs[datapairs['charge_transfer'] != None]



# # Sorts the frames byt their name so that they are plotted alphabetically
# datapairs = datapairs.sort_values(["ion_name", "sort"])

# # Sets the dataframe column type to float
# datapairs.electrostatics = datapairs.electrostatics.astype(float)
# datapairs.exchange = datapairs.exchange.astype(float)
# datapairs.induction = datapairs.induction.astype(float)
# datapairs.dispersion = datapairs.dispersion.astype(float)
# datapairs.total = datapairs.total.astype(float)
# datapairs.charge_transfer = datapairs.charge_transfer.astype(float)
# datapairs.dhomo = datapairs.dhomo.astype(float)
# datapairs.dlumo = datapairs.dlumo.astype(float)
# datapairs.bhomo = datapairs.bhomo.astype(float)
# datapairs.blumo = datapairs.blumo.astype(float)
# datapairs.plotsize = datapairs.plotsize.astype(int)

# datapairs = pd.merge(datapairs, homodatapairs, on=['base_name', 'ion_name'])

# Eh2eV(datapairs, ["lumo", "homo", "dhomo", "dlumo", "bhomo", "blumo"])
# Eh2eV(lonepairs, ["lumo", "homo"])

# warnings.filterwarnings("ignore")
# for i in range(len(datapairs.base_name)):
#     for j in range(len(lonepairs.homo)):

#         if datapairs.base_name[i] == lonepairs.base_name[j]:
#             datapairs.dhomo[i] = datapairs.homo[i] - lonepairs.homo[j]
#             datapairs.dlumo[i] = datapairs.lumo[i] - lonepairs.lumo[j]
#             datapairs.bhomo[i] = lonepairs.homo[j]
#             datapairs.blumo[i] = lonepairs.lumo[j]

# for i in range(len(datapairs.base_name)):
#     if math.isnan(datapairs.total[i]) == False:
#         datapairs.plotsize[i] = abs(int(datapairs.total[i]))

Export data to csv for future easy import

# data.to_csv(r'./IonsData/pd.data.csv')
# datapairs.to_csv(r'./IonsData/pd.datapairs.csv')

(new) Import and process data

data = pd.read_csv("./IonsData/pd.data.csv")
datapairs = pd.read_csv("./IonsData/pd.datapairs.csv")

data = data.sort_values(["ion_name", "sort"])

datapairs.electrostatics = datapairs.electrostatics.astype(float)
datapairs.exchange = datapairs.exchange.astype(float)
datapairs.induction = datapairs.induction.astype(float)
datapairs.dispersion = datapairs.dispersion.astype(float)
datapairs.total = datapairs.total.astype(float)
datapairs.charge_transfer = datapairs.charge_transfer.astype(float)
datapairs.dhomo = datapairs.dhomo.astype(float)
datapairs.dlumo = datapairs.dlumo.astype(float)
datapairs.bhomo = datapairs.bhomo.astype(float)
datapairs.blumo = datapairs.blumo.astype(float)
datapairs.plotsize = datapairs.plotsize.astype(int)

data.electrostatics = data.electrostatics.astype(float)
data.exchange = data.exchange.astype(float)
data.induction = data.induction.astype(float)
data.dispersion = data.dispersion.astype(float)
data.total = data.total.astype(float)
data.charge_transfer = data.charge_transfer.astype(float)
data.dhomo = data.dhomo.astype(float)
data.dlumo = data.dlumo.astype(float)
data.bhomo = data.bhomo.astype(float)
data.blumo = data.blumo.astype(float)
data.plotsize = data.plotsize.astype(int)

Lone Ion Plots

HOMO and $\Delta$HOMO

Please disregard Benzene-Ni

def plothomo(df, ionclass):
    if ionclass == "metal":
        widthval = 500
    else:
        widthval = 1000
    fig = px.bar(df[df.ion_class == ionclass],
                 x="base_name",
                 y="homo",
                 facet_col='ion_name',
                 labels={'homo':'Energy (eV)', 'base_name':""},
                 height=400,
                 width=widthval,
                 barmode="overlay",
                 opacity=1,
                )
    for i in range(len(fig.data)):
        fig.data[i].name="HOMO"
    fig2 = px.bar(df[df.ion_class == ionclass],
                 x="base_name",
                 y="dhomo",
                 facet_col='ion_name',
                 labels={'homo':'Energy (eV)', 'base_name':""},
                 height=400,
                 width=widthval,
                 barmode="overlay",
                 opacity=1,
                 color_discrete_map={"dhomo": "rgb(255,119,119)"}
                )

    for i in range(len(fig2.data)):
        fig2.data[i].name="ΔHOMO"
        
    fig.add_trace(fig2.data[0])
    fig.add_trace(fig2.data[1])
    if ionclass != "metal":
        fig.add_trace(fig2.data[2])
        fig.add_trace(fig2.data[3])
        fig.add_trace(fig2.data[4])
        fig.add_trace(fig2.data[5])
        if ionclass != "cation":
            fig.add_trace(fig2.data[6])
    fig.data[1].showlegend = True
    fig.data[-1].showlegend = True
    fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
    return(fig.show())


for ionclass in classnames:
    plothomo(data, ionclass)
    

SAPT Energies

Excluding Benzene-Ni

def saptdf(df):
    newdf = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'energy', 'energytype'])
    for i in list(range(len(df))):
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.total[i],"energytype" : "total"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.electrostatics[i],"energytype" : "electrostatics"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.exchange[i],"energytype" : "exchange"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.induction[i],"energytype" : "induction"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.dispersion[i],"energytype" : "dispersion"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.charge_transfer[i],"energytype" : "charge_transfer"},ignore_index=True)
        
    newdf.sort_values(["ion_name", "sort"])
    return(newdf)

def plotsapt(df, ionclass):
    dfint = saptdf(df)
    if ionclass == "metal":
        widthval = 550
    else:
        widthval = 1000
        
    fig = px.scatter(dfint[dfint.ion_class == ionclass],
                 x="base_name",
                 y="energy",
                 facet_col='ion_name',
                 labels={'energy':'Energy (KJ/mol)', 'base_name':""},
                 height=400,
                 width=widthval,
#                  barmode="overlay",
#                  opacity=0.5,
                 color="energytype",
                )
    fig.layout.showlegend = True
    fig.update_traces(opacity=0.7, marker=dict(size=10, line=dict(color='black', width=1)))
    fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
    fig.for_each_trace(lambda t: t.update(name=t.name.replace("energytype=","")))
    return(fig.show())


for ionclass in classnames:
    plotsapt(data, ionclass)

   trendlist = ["ion_name", "ion_class", "base_name"]

Correlation Plot

Coloured by Ion Name

Excluding Pd$^{2+}$, Ni$^{2+}$ and hard ions(Cl$^-$, Na$^+$ and Li$^+$)

dplusel is dispersion + electrostatics

# Feature Plot


datacorr = data
datacorr['dplusel'] = np.add(datacorr['electrostatics'], datacorr['dispersion'])

# Correlation plot
fig = px.scatter_matrix(datacorr[datacorr["ion_class"] != "metal"][datacorr["ion_name"] != "Na"][datacorr["ion_name"] != "Li"][datacorr["ion_name"] != "Cl"],
                    dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
                    color="ion_name",
                    hover_data=['ion_name', 'base_name'],
                    width=1000,
                    height=1000)

fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()

Coloured by Ion Class (Anion/Cation)

Excluding Pd$^{2+}$, Ni$^{2+}$ and hard ions(Cl$^-$, Na$^+$ and Li$^+$)

dplusel is dispersion + electrostatics

datacorr = data
datacorr['dplusel'] = np.add(datacorr['electrostatics'], datacorr['dispersion'])
datacorr['dpluselplusex'] = np.add(datacorr['dplusel'], datacorr['exchange'])

# Correlation plot
fig = px.scatter_matrix(datacorr[datacorr["ion_class"] != "metal"][datacorr["ion_name"] != "Na"][datacorr["ion_name"] != "Li"][datacorr["ion_name"] != "Cl"],
                    dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
                    color="ion_class",
                    hover_data=['ion_name', 'base_name'],
                    width=1000,
                    height=1000)

fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()

Ion Pair Plots

HOMO and $\Delta$HOMO

def plothomopairs(df):
    fig = px.bar(df,
                 x="base_name",
                 y="homo",
                 facet_col='ion_name',
                 labels={'homo':'Energy (eV)', 'base_name':""},
                 width=600,
                 height=400,
                 barmode="overlay",
                 opacity=1,
                )
    for i in range(len(fig.data)):
        fig.data[i].name="HOMO"
    fig2 = px.bar(df,
                 x="base_name",
                 y="dhomo",
                 facet_col='ion_name',
                 labels={'homo':'Energy (eV)', 'base_name':""},
                 height=400,
                 barmode="overlay",
                 opacity=1,
                 color_discrete_map={"dhomo": "rgb(31,120,180)"}
                )
    for i in range(len(fig2.data)):
        fig2.data[i].name="ΔHOMO"
    fig.add_trace(fig2.data[0])
    fig.add_trace(fig2.data[1])

    fig.data[1].showlegend = True
    fig.data[-1].showlegend = True
    fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
    return(fig.show())


plothomopairs(datapairs)
        
        

SAPT Energy

def saptdfpairs(df):
    newdf = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'energy', 'energytype'])
    for i in list(range(len(df))):
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.total[i],"energytype" : "total"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.electrostatics[i],"energytype" : "electrostatics"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.exchange[i],"energytype" : "exchange"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.induction[i],"energytype" : "induction"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.dispersion[i],"energytype" : "dispersion"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i], "ion_name" : df.ion_name[i],"energy" : df.charge_transfer[i],"energytype" : "charge_transfer"},ignore_index=True)
        
    newdf.sort_values(["ion_name", "sort"])
    return(newdf)

def plotsaptpairs(df):
    dfint = saptdfpairs(df)
    if ionclass == "metal":
        widthval = 550
    else:
        widthval = 600
        
    fig = px.scatter(dfint,
                 x="base_name",
                 y="energy",
                 facet_col='ion_name',
                 labels={'homo':'Energy (eV)', 'base_name':""},
                 width=600,
                 height=400,
#                  width=widthval,
#                  barmode="overlay",
                 color="energytype",
                )
    fig.layout.showlegend = True
    fig.update_traces(opacity=0.7, marker=dict(size=12, line=dict(color='black', width=1)))
    fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
    fig.for_each_trace(lambda t: t.update(name=t.name.replace("energytype=","")))
    return(fig.show())


plotsaptpairs(datapairs)

Correlation Plot

Both

Correlation Plot

Ion pairs and Individual ions combined

Excluding Pd$^{2+}$, Ni$^{2+}$, and hard ions(Cl$^-$, Na$^+$ and Li$^+$)

databoth = pd.DataFrame(columns=["base_name", "ion_name", "electrostatics", "exchange", "induction", "dispersion", "total", "charge_transfer", "dhomo", "sorp"])
for i in range(len(list(data.base_name))):
    newdata = {"base_name":data.base_name[i], "ion_name":data.ion_name[i], "electrostatics":data.electrostatics[i], "exchange":data.exchange[i], "induction":data.induction[i], "dispersion":data.dispersion[i], "total":data.total[i], "charge_transfer":data.charge_transfer[i], "dhomo":data.dhomo[i], "sorp":"s"}
    databoth = databoth.append(newdata, ignore_index=True)
for i in range(len(list(datapairs.base_name))):
    newdata = {"base_name":datapairs.base_name[i], "ion_name":datapairs.ion_name[i], "electrostatics":datapairs.electrostatics[i], "exchange":datapairs.exchange[i], "induction":datapairs.induction[i], "dispersion":datapairs.dispersion[i], "total":datapairs.total[i], "charge_transfer":datapairs.charge_transfer[i], "dhomo":datapairs.dhomo[i], "sorp":"p"}
    databoth = databoth.append(newdata, ignore_index=True)
    
    
databoth = databoth[databoth["ion_name"] != "Ni"][databoth["ion_name"] != "Pd"][databoth["ion_name"] != "Na"][databoth["ion_name"] != "Cl"][databoth["ion_name"] != "Li"][databoth["ion_name"] != "Na"]
databoth['dplusel'] = np.add(databoth['electrostatics'], databoth['dispersion'])
databoth['dpluselplusex'] = np.add(databoth['dplusel'], databoth['exchange'])


fig = px.scatter_matrix(databoth,
                    dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
                    color="sorp",
                    hover_data=['ion_name', 'base_name'],
                    width=1000,
                    height=1000)

fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()

Katya's Suggestions

1) Correlations

Except Li (Na is missing for some reason) the electrostatics follows the total interaction energy. Could you please extract the gradient from the electrostatics-total graph. I believe that Induction + dispersion + exchange might be summing to within a few kJ mol-1.... Ultimately, electrostatics is the force that drives this interaction... Li definitely makes the bond covalent and hence such a huge difference between electrostatics and total energy.

Total - Electrostatics

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$

data['dplusel'] = np.add(data['electrostatics'], data['dispersion'])
data['dpluselplusex'] = np.add(data['dplusel'], data['exchange'])

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
                 y="total",
                 x="electrostatics",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("total = " + str(results.px_fit_results.iloc[0].params[1]) +"* Electrostatics " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.9030553646414368
total = 0.9856406966724351* Electrostatics -6.433560778791372

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$, Na$^+$, Cl$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
                 y="total",
                 x="electrostatics",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Total = " + str(results.px_fit_results.iloc[0].params[1]) +" * Electrostatics + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.8904670996931319
Total = 0.8736926799703468 * Electrostatics + -7.35623462676404

Total - Dispersion + Electrostatics

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
                 y="total",
                 x="dplusel",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("total = " + str(results.px_fit_results.iloc[0].params[1]) +"* (Electrostatics + Dispersion) " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.6289566526512076
total = 0.743313858961744* (Electrostatics + Dispersion) 3.291513720534722

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$, Na$^+$, Cl$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
                 y="total",
                 x="dplusel",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("total = " + str(results.px_fit_results.iloc[0].params[1]) +"* (Electrostatics + Dispersion) " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.7681970967915128
total = 0.6042561976435759* (Electrostatics + Dispersion) 1.9471283622225952

Total - Dispersion + Electrostatics + Exchange

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
                 y="total",
                 x="dpluselplusex",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Total = " + str(results.px_fit_results.iloc[0].params[1]) +" * (Electrostatics + Dispersion + Exchange) + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.7918257575866938
Total = 1.7648848320532988 * (Electrostatics + Dispersion + Exchange) + -5.364755450287623

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$, Na$^+$, Cl$^+$

fig = px.scatter(data[datacorr.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
                 y="total",
                 x="dpluselplusex",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Total = " + str(results.px_fit_results.iloc[0].params[1]) +" * (Electrostatics + Dispersion + Exchange) + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.9424319306598973
Total = 1.4124998291175728 * (Electrostatics + Dispersion + Exchange) + -6.249230648555187

2) Further Correlations

Electrostatics and exchange also correlate with each other, and so do a) dispersion and exchange and b) induction and electrostatics. Extracting their gradients would be good too. Surprisingly, dispersion correlates with exchange without the need to add electrostatics! Although the latter does improve the correlation.... The dispersion-exchange correlation could be compared to Anh's data.... I think this is the reason why her DFT works so well...

Electrostatics - Exchange

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
                 y="electrostatics",
                 x="exchange",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Electrostatics = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.5118609630731161
Electrostatics = -1.063969275341984 * Exchange + 9.512563076218635

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$, Na$^+$, Cl$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
                 y="electrostatics",
                 x="exchange",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Electrostatics = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.5316140446902322
Electrostatics = -0.8780190254123053 * Exchange + 8.361633853052801

Dispersion - Exchange

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$

fig = px.scatter(datacorr[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
                 y="dispersion",
                 x="exchange",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Dispersion = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.37000294355481467
Dispersion = -0.4877974627087508 * Exchange + -6.003526161695268

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$, Na$^+$, Cl$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
                 y="dispersion",
                 x="exchange",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Dispersion = " + str(results.px_fit_results.iloc[0].params[1]) +" * Exchange + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.7464544651613487
Dispersion = -0.6360920365032324 * Exchange + -4.788495179125986

Induction - Electrostatics

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"],
                 y="induction",
                 x="electrostatics",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Induction = " + str(results.px_fit_results.iloc[0].params[1]) +" * Electrostatics + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.6506407144336758
Induction = 0.4998111176981335 * Electrostatics + -4.337012515738229

Excluding Ni$^{2+}$, Pd$^{2+}$, Li$^+$, Na$^+$, Cl$^+$

fig = px.scatter(data[data.ion_name != "Ni"][data.ion_name != "Pd"][data.ion_name != "Li"][data.ion_name != "Na"][data.ion_name != "Cl"],
                 y="induction",
                 x="electrostatics",
#                  color="sorp",
                 height=400,
                 width=600,
                 hover_data=['ion_name', 'base_name'],
                 trendline="ols",
#                      size='plotsize'
                 )
fig.update_layout(legend_orientation="h")
fig.show()
results = px.get_trendline_results(fig)
print("R^2 = " + str(results.px_fit_results.iloc[0].rsquared))
print("Induction = " + str(results.px_fit_results.iloc[0].params[1]) +" * Electrostatics + " + str(results.px_fit_results.iloc[0].params[0]))

R^2 = 0.5906389608317781
Induction = 0.26427779326417417 * Electrostatics + -6.42340309311744

3) Correlation of ONLY ionpairs and their corresponding lone counterparts

I would like to know where the ion pairs are located on 3.1.1 plots. Could you please include only those of the lowest energy? Would you be able to make SAPT energy graphs combining the ion pairs and the results for individual ions included in these ion pairs?

databoth2 = pd.DataFrame(columns=["base_name", "ion_name", "electrostatics", "exchange", "induction", "dispersion", "total", "charge_transfer", "dhomo", "sorp"])
for i in range(len(list(data.base_name))):
    newdata = {"base_name":data.base_name[i], "ion_name":data.ion_name[i], "electrostatics":data.electrostatics[i], "exchange":data.exchange[i], "induction":data.induction[i], "dispersion":data.dispersion[i], "total":data.total[i], "charge_transfer":data.charge_transfer[i], "dhomo":data.dhomo[i], "sorp":"s"}
    databoth2 = databoth2.append(newdata, ignore_index=True)
for i in range(len(list(datapairs.base_name))):
    newdata = {"base_name":datapairs.base_name[i], "ion_name":datapairs.ion_name[i], "electrostatics":datapairs.electrostatics[i], "exchange":datapairs.exchange[i], "induction":datapairs.induction[i], "dispersion":datapairs.dispersion[i], "total":datapairs.total[i], "charge_transfer":datapairs.charge_transfer[i], "dhomo":datapairs.dhomo[i], "sorp":"p"}
    databoth2 = databoth2.append(newdata, ignore_index=True)
    
databoth2 = databoth2[databoth2.ion_name != "BF2(CF3)2"][databoth2.ion_name != "C(CN)3"][databoth2.ion_name != "N(SOOCF3)2"][databoth2.ion_name != "N(NC)2"][databoth2.ion_name != "C1MPyr"][databoth2.ion_name != "NMe4"][databoth2.ion_name != "PMe4"][databoth2.ion_name != "Na"][databoth2.ion_name != "Pd"][databoth2.ion_name != "Ni"][databoth2.ion_name != "tempo"]
databoth2['dplusel'] = np.add(databoth2['electrostatics'], databoth2['dispersion'])
databoth2['dpluselplusex'] = np.add(databoth2['dplusel'], databoth2['exchange'])

databoth2

fig = px.scatter_matrix(databoth2,
                    dimensions=["total", "electrostatics", "exchange", "induction", "dispersion", "charge_transfer", "dhomo", "dplusel"],
                    color="sorp",
                    hover_data=['ion_name', 'base_name'],
                    width=1000,
                    height=1000)

fig.update_traces(diagonal_visible=False, showupperhalf=True)
fig.update_layout(legend_orientation="h")
fig.show()

6) Induction - Electrostatics SAPT comparison

Ni and Pd form covalent bonds by the look of it. In both cases it is the induction that dominates the interaction (hence, the symmetry is crucial). Electrostatics is stronger in Ni but is compensated by high exchange... So induction is the key whereas in ILs pure electrostatics is the key! It must come down to the orbital symmetry I guess. This is very cool though as clearly transition metals favour certain sites whereas the problem of ions is that they lead to non-directional interaction... So, the next goal would be to identify ions that have larger induction over electrostatics (if it is even possible).

This problem seems to be dependent on the molecule with which the ion is interacting and unfortunately we can only really say that BF4$^+$ and Li$^-$ have the highest number of incidence where induction is greater than electrostatics.

def saptdf(df):
    newdf = pd.DataFrame(columns=['sort', 'base_name', 'ion_class', 'ion_name', 'energy', 'energytype'])
    for i in list(range(len(df))):
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.electrostatics[i],"energytype" : "electrostatics"},ignore_index=True)
        newdf = newdf.append({"sort": df.sort[i],"base_name" : df.base_name[i],"ion_class" : df.ion_class[i],"ion_name" : df.ion_name[i],"energy" : df.induction[i],"energytype" : "induction"},ignore_index=True)
  
    newdf.sort_values(["ion_name", "sort"])
    return(newdf)

def plotsapt(df, ionclass):
    dfint = saptdf(df)
    if ionclass == "metal":
        widthval = 550
    else:
        widthval = 1000
        
    fig = px.scatter(dfint[dfint.ion_class == ionclass],
                 x="base_name",
                 y="energy",
                 facet_col='ion_name',
                 labels={'energy':'Energy (KJ/mol)', 'base_name':""},
                 height=400,
                 width=widthval,
#                  barmode="overlay",
#                  opacity=0.5,
                 color="energytype",
                )
    fig.layout.showlegend = True
    fig.update_traces(opacity=0.7, marker=dict(size=10, line=dict(color='black', width=1)))
    fig.for_each_annotation(lambda a: a.update(text=a.text.replace("ion_name=", "")))
    fig.for_each_trace(lambda t: t.update(name=t.name.replace("energytype=","")))
    return(fig.show())


for ionclass in classnames:
    plotsapt(data, ionclass)